CITY_US
I Предварительный анализ данных
if (interactive() && !str_ends(getwd(), "R/Statistical_Analysis_2022/CITY_US/Tkachenko")) {
setwd("R/Statistical_Analysis_2022/CITY_US/Tkachenko")
}
data <- read_excel("../CITY_US.STD/CITY_shortname.xls")
data[data == "NA"] <- NA
data[, -(1:2)] <- data.frame(lapply(data[, -(1:2)], as.numeric))
fullnames <- names(read_excel("../CITY_US.STD/CITY.xls"))1 Разобраться в том, что означают признаки.
print(fullnames)## [1] "CITY City"
## [2] "STATE State"
## [3] "AREA Land area, 1990 (sq.miles)"
## [4] "R_AREA Land area, 1990 (sq.miles), ranks"
## [5] "POP92 Population,1992"
## [6] "R_POP92 Population,1992, ranks"
## [7] "POP80 Population,1980"
## [8] "R_POP80 Population,1980, ranks"
## [9] "POP_CH Population growth rate,1980-1992"
## [10] "R_POP_CH Population growth rate,1980-1992, ranks"
## [11] "POPDEN Population per square mile, 1992"
## [12] "R_POPDEN Population per square mile, 1992, ranks"
## [13] "OLD % older 65 years"
## [14] "R_OLD % older 65 years, ranks"
## [15] "BLACK Black population,1990"
## [16] "R_BLACK Black population,1990, ranks"
## [17] "BLACK% % Black population,1990"
## [18] "R_BLACK% % Black population,1990, ranks"
## [19] "ASIAN Asian or Pacific Islander population,1990"
## [20] "R_ASIAN Asian or Pacific Islander population,1990, ranks"
## [21] "ASIAN% % Asian or Pacific Islander population,1990"
## [22] "R_ASIAN% % Asian or Pacific Islander population,1990, ranks"
## [23] "HISP Hispanic population, 1990"
## [24] "R_HISP Hispanic population, 1990, ranks"
## [25] "HISP% % Hispanic population, 1990"
## [26] "R_HISP% % Hispanic population, 1990, ranks"
## [27] "BORN_F % persons of foreign born"
## [28] "R_BORN_F % persons of foreign born, ranks"
## [29] "LANG % of persons, speaking language other than English at home, 1990"
## [30] "R_LANG % of persons, speaking language other than English at home, 1990, ranks"
## [31] "HH1 % 1-person households, 1990"
## [32] "R_HH1 % 1-person households, 1990, ranks"
## [33] "FAMIL1 % one-parent family households, 1990"
## [34] "R_FAMIL1 % one-parent family households, 1990, ranks"
## [35] "DEATH Infant death rate per 1000 live births,1988"
## [36] "R_DEATH Infant death rate per 1000 live births,1988, ranks"
## [37] "CRIME Crime rate per 100000 population, 1991"
## [38] "R_CRIME Crime rate per 100000 population, 1991, rate"
## [39] "SCHOOL % of elementary and high school enrollment in public schools, 1990"
## [40] "R_SCHOOL % of elementary and high school enrollment in public schools, 1990, ranks"
## [41] "DEGREE % of adults with a bachelor's degree or higher, 1990"
## [42] "R_DEGREE % of adults with a bachelor's degree or higher, 1990, ranks"
## [43] "INCOME Median household income, 1989"
## [44] "R_INCOME Median household income, 1989, ranks"
## [45] "ASSIST % of households, receiving public assistance,1989"
## [46] "R_ASSIST % of households, receiving public assistance,1989, ranks"
## [47] "POVERT % of persons below poverty level, 1989"
## [48] "R_POVERT % of persons below poverty level, 1989, ranks"
## [49] "OLB_BIL % of housing units built 1939 or earlier, 1990"
## [50] "R_OLD_BI % of housing units built 1939 or earlier, 1990, ranks"
## [51] "OWNER Median value of specified owner-occupied housing units ($), 1990"
## [52] "R_OWNER Median value of specified owner-occupied housing units ($), 1990, ranks"
## [53] "RENTER % renter-occupied housing units,1990"
## [54] "R_RENTER % renter-occupied housing units,1990, ranks"
## [55] "GROSS Median gross rent of specified renter-occupied housing units ($), 1990"
## [56] "R_GROSS Median gross rent of specified renter-occupied housing units ($), 1990, ranks"
## [57] "CONDOM % condominimum occupied housing units, 1990"
## [58] "R_CONDOM % condominimum occupied housing units, 1990, ranks"
## [59] "TRANSP % of workers using public transportation"
## [60] "R_TRANSP % of workers using public transportation, ranks"
## [61] "UNEMP Unemployment rate, 1991"
## [62] "R_UNEMP Unemployment rate, 1991, ranks"
## [63] "LABOR Labor force - % change 1980-1990"
## [64] "R_LABOR Labor force - % change 1980-1990, ranks"
## [65] "LAB_F Female civilian labor force participation rate, 1990"
## [66] "R_LAB_F Female civilian labor force participation rate, 1990, ranks"
## [67] "MANLAB % of civilian labor force emploied in manufacturing, 1990"
## [68] "R_MANLAB % of civilian labor force emploied in manufacturing, 1990, ranks"
## [69] "TAXE City government taxes per capita, 1991"
## [70] "R_TAXE City government taxes per capita, 1991, ranks"
## [71] "TEMPER Average daily temperature in July"
## [72] "R_TEMPER Average daily temperature in July, ranks"
## [73] "PRECEP Annual precipitation"
## [74] "R_PRECEP Annual precipitation"2 Отобрать признаки
names_interesting <- c("AREA", "POP80", "POP92", "POPDEN", "CRIME", "BORN_F", "POVERT", "INCOME", "UNEMP", "TEMPER")
data <- data %>% select(all_of(c("CITY", "STATE", names_interesting)))
print(head(data))## # A tibble: 6 × 12
## CITY STATE AREA POP80 POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 NEW_… NY 309. 7.07e6 7.31e6 23671 9236 28.4 19.3 29823 8.6 76.8
## 2 LOS_… CA 469. 2.97e6 3.49e6 7436 9730 38.4 18.9 30925 9 74.3
## 3 CHIC… IL 227. 3.01e6 2.77e6 12185 NA 16.9 21.6 26301 8.4 75.1
## 4 HOUS… TX 540. 1.60e6 1.69e6 3131 10824 17.8 20.7 26261 6.1 83.5
## 5 PHIL… PA 135. 1.69e6 1.55e6 11492 6835 6.6 20.3 24603 8 76.7
## 6 SAN_… CA 324 8.76e5 1.15e6 3546 8537 20.9 13.4 10 6.2 713 Определить вид признаков
Город и штат качественные, остальные количественные, ранги были порядковыми.
find_mode_freq <- function(x) {
x <- x[!is.na(x)]
return(max(tabulate(match(x, x))))
}
print(data %>% summarise(across(
all_of(names_interesting),
find_mode_freq
)))## # A tibble: 1 × 10
## AREA POP80 POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
## <int> <int> <int> <int> <int> <int> <int> <int> <int> <int>
## 1 1 1 1 2 1 4 2 1 5 3print(sort(data$UNEMP))## [1] 2.3 3.2 3.8 3.9 4.1 4.4 4.6 4.7 4.7 4.8 4.8 5.0 5.0 5.0 5.0
## [16] 5.1 5.1 5.1 5.3 5.4 5.4 5.4 5.4 5.4 5.6 5.7 5.9 5.9 5.9 6.0
## [31] 6.0 6.1 6.1 6.1 6.1 6.2 6.2 6.4 6.4 6.5 6.7 6.7 6.7 6.8 6.9
## [46] 6.9 7.0 7.1 7.2 7.2 7.3 7.3 7.3 7.5 7.6 7.7 7.7 7.7 8.0 8.1
## [61] 8.4 8.4 8.5 8.6 9.0 9.0 9.4 9.5 9.6 10.0 10.4 10.6 10.7 11.0 11.9
## [76] 12.6 13.1Все количественные буду считать непрерывными. возможно UNEMP непрерывный с плохой точностью.
4 не актуально
5 Построить matrix plot
if (interactive()) pdf("ggpairs_unedited.pdf")
ggpairs(
data[, -(1:2)],
lower = list(continuous = wrap("points", alpha = 0.5, size = 0.3)),
diag = list(continuous = "barDiag")
)
if (interactive()) dev.off()7 outliers
Убираю outliers:
Помечаю некорректные данные в INCOME как NA. Удаляю город из Аляски за плотность населения. Флорида выделсется на BORN_F-INCOME Гаваи выделяются низким уровнем безработицы. странно?
data$INCOME[data$INCOME < 100] <- NA
data <- data %>% filter(STATE != "AK")6 Несимметричные распределения
Функция, которая логарифмирует, если это сделает выборку симметричнее
log_asymmetric <- function(x) {
if (skewness(x, na.rm = TRUE) < abs(skewness(log(x), na.rm = TRUE))) {
print("default")
return(x)
} else {
print("logged")
return(log(x))
}
}Автоматически логарифмирую то что имеет асимметрию и длинный хвост справа
data_logged <- data %>%
mutate(across(all_of(names_interesting), log_asymmetric))## [1] "logged"
## [1] "logged"
## [1] "logged"
## [1] "logged"
## [1] "logged"
## [1] "logged"
## [1] "default"
## [1] "default"
## [1] "logged"
## [1] "default"if (interactive()) pdf("ggpairs_logged.pdf")
ggpairs(
data_logged[, -(1:2)],
lower = list(continuous = wrap("points", alpha = 0.5, size = 0.3)),
diag = list(continuous = "barDiag")
)
if (interactive()) dev.off()8 однородность
выглядит однородно.
9 не актуально
10 всякие характеристики
print_characteristics <- function(x) {
list(
mean = mean(x, na.rm = TRUE),
var = var(x, na.rm = TRUE),
skewness = skewness(x, na.rm = TRUE),
kurtosis = kurtosis(x, na.rm = TRUE) - 3
)
}
sapply(data_logged %>% select(all_of
(names_interesting)), print_characteristics)## AREA POP80 POP92 POPDEN CRIME BORN_F
## mean 4.754638 12.9069 13.01226 8.257586 9.206558 1.959747
## var 0.6658984 0.5142076 0.4464288 0.5099096 0.06798296 0.8435754
## skewness 0.1317646 1.453879 1.743773 0.1015321 0.1476689 0.3081197
## kurtosis -0.3907363 3.033026 3.885752 -0.1531417 0.08872069 -0.7431516
## POVERT INCOME UNEMP TEMPER
## mean 18.11842 25282.48 1.873674 77.60132
## var 34.68872 14321837 0.09967048 37.59853
## skewness 0.2243883 -0.2196584 -0.2186434 -0.1426746
## kurtosis -0.3539569 -0.6467894 0.6594639 0.7472063II О виде распределений и о сравнении распределений
1 выполнить первое второе задание в логичном порядке
2 анализ вида распределения признаков
2. Описать распределения признаков визуально (с помощью нормальной бумаги и PP-plot) и по критерию хи-квадрат и другим критериям.
Сюда входит: normal probability plot (что это такое?), проверка по критериям Лиллиефорса, AD, хи-квадрат, Шапиро-Уилка. По критерию хи-квадрат, а также визуально по PP-plot можно проверить и гипотезы о согласии с другими распределениями, например, логнормальным.
Рассматриваем нормальность
все плотности.
if (interactive()) pdf(file = "all_density_plots.pdf")
ggplot(melt(data_logged, id = c("CITY", "STATE")), aes(value)) +
theme_bw() +
geom_histogram(aes(y = ..density..), bins = 15) +
# geom_density() +
facet_wrap(~variable, scales = "free") +
geom_line(aes(y = dnorm(value,
mean = tapply(value, variable, mean, na.rm = TRUE)[PANEL],
sd = tapply(value, variable, sd, na.rm = TRUE)[PANEL]
)), color = "red") +
theme(legend.position = "bottom") +
labs(x = "", y = "")
if (interactive()) dev.off()Normal probability plot
if (interactive()) pdf(file = "normal_probability_plot.pdf")
ggplot(melt(data_logged, id = c("CITY", "STATE")), aes(sample = value)) +
stat_qq_point(size = 2) +
geom_abline() +
facet_wrap(~variable, scales = "free")
if (interactive()) dev.off()PP-plot
if (interactive()) pdf(file = "PP_plot.pdf")
ggplot(melt(data_logged, id = c("CITY", "STATE")), aes(sample = value)) +
stat_pp_point(size = 2) +
facet_wrap(~variable, scales = "free") +
stat_pp_line()
if (interactive()) dev.off()Всякие тесты на нормальность
data_logged$F_TEMPER <- NULL
test_all <- function(data, test, ...) {
print(test(1:10, ...)$method)
print(data %>% summarise(across(
everything(),
function(x) {
test(x, ...)$p.value
}
)))
}
test_all(data_logged[, -(1:2)], shapiro.test)## [1] "Shapiro-Wilk normality test"
## # A tibble: 1 × 10
## AREA POP80 POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.757 0.0000113 0.000000310 0.673 0.725 0.0892 0.675 0.463 0.487 0.368test_all(data_logged[, -(1:2)], lillie.test)## [1] "Lilliefors (Kolmogorov-Smirnov) normality test"
## # A tibble: 1 × 10
## AREA POP80 POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.536 0.0118 0.000633 0.291 0.580 0.148 0.646 0.866 0.746 0.766test_all(data_logged[, -(1:2)], ad.test)## [1] "Anderson-Darling normality test"
## # A tibble: 1 × 10
## AREA POP80 POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.492 0.000167 0.000000752 0.370 0.620 0.0362 0.830 0.692 0.697 0.371test_all(data_logged[, -(1:2)], pearson.test)## [1] "Pearson chi-square normality test"
## # A tibble: 1 × 10
## AREA POP80 POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.697 0.0377 0.0179 0.664 0.728 0.265 0.729 0.857 0.897 0.0626Критерии отвергают нормальность распределения населения. В остальных признаках отклонения от нормального (логнормального для ранее логарифмированных признаков) распределения тестами не обнаружено.
1. Описать разницу между численностью населения в 1980 и 1992 годах визуально (с помощью ящиков с усами) и по критериям. Добавить признак, который разбивает города на “холодные” и “жаркие”. Сравнить их по характеристикам.
добавим признак
temper_split <- median(data_logged$TEMPER)
data_logged$F_TEMPER <- factor(ifelse(data_logged$TEMPER < temper_split, "cold", "hot"))3 Сначала имеет смысл посмотреть на сравнение распределений в группах с помощью ящиков с усами
if (interactive()) pdf("boxplot_pop.pdf")
ggplot(
melt(data_logged, id.vars = "F_TEMPER", measure.vars = c("POP80", "POP92")),
aes(x = variable, y = value, color = F_TEMPER)
) +
geom_boxplot() +
scale_color_manual(values = c("blue", "red"))
if (interactive()) dev.off()
if (interactive()) pdf("boxplot_pop_all.pdf")
ggplot(
melt(data_logged, id.vars = "F_TEMPER", measure.vars = c("AREA", "POP80", "POPDEN", "CRIME", "BORN_F", "POVERT", "INCOME", "UNEMP", "TEMPER")),
aes(x = variable, y = value, color = F_TEMPER)
) +
geom_boxplot() +
facet_wrap(~variable, scales = "free") +
scale_color_manual(values = c("blue", "red"))## Warning: Removed 11 rows containing non-finite values (`stat_boxplot()`).
if (interactive()) dev.off()Сравнение признаков при разбиении городов на группы по температуре
4 t-test
Разделение на группы было сделано по медиане, поэтому дизайн сбалансированный, и для t-теста можно не выяснять равенство дисперсий.
Но из интереса применю тест Фишера, а перед этим проверю группы на нормальность распределения.
Признак CRIME
Из ящиков с усами ожидаю отвержение нулевых гипотез в тестах.
hot <- (data_logged %>%
filter(F_TEMPER == "hot"))$CRIME
cold <- (data_logged %>%
filter(F_TEMPER == "cold"))$CRIMEcat(shapiro.test(hot)$p.value, shapiro.test(cold)$p.value)## 0.5012446 0.954743cat(pearson.test(hot)$p.value, pearson.test(cold)$p.value)## 0.1314607 0.4171623Тесты на нормальность не обнаружили ненормальности в распределениях. Применю тест Фишера
var.test(hot, cold)##
## F test to compare two variances
##
## data: hot and cold
## F = 1.6011, num df = 37, denom df = 36, p-value = 0.1608
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.8269838 3.0904136
## sample estimates:
## ratio of variances
## 1.601064Тест Фишера не отверг гипотезу о равенстве дисперсий, но это не так важно, потому что в t-test дизайн сбалансирован и формула не изменится.
t.test(cold, hot)##
## Welch Two Sample t-test
##
## data: cold and hot
## t = -3.2207, df = 70.063, p-value = 0.00194
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.29506315 -0.06938105
## sample estimates:
## mean of x mean of y
## 9.114233 9.296455t-test уверенно отвергает гипотезу о равенстве матожиданий, как и ожидалось.
5 непараметрический тест Вилкоксона.
Если бы дизайн был несбалансирован или не было бы нормальности распределений и нас не удовлетворяла асимптотичность t-критерия или нас интересовала другая характеристика положение, то мы бы применили критерий Вилкоксона.
wilcox.test(cold, hot, paired = FALSE)##
## Wilcoxon rank sum exact test
##
## data: cold and hot
## W = 413, p-value = 0.001867
## alternative hypothesis: true location shift is not equal to 0Он тоже отвергает схожесть распределений.
7 критерий Колмогорова-Смирнова
ks.test(cold, hot)##
## Exact two-sample Kolmogorov-Smirnov test
##
## data: cold and hot
## D = 0.38762, p-value = 0.004648
## alternative hypothesis: two-sidedТест Колмогорова-Смирнова обнаруживает разницу в форме распределений.
Признак INCOME
для этого признака многое повторится, но я ожидаю неотвержения нулевых гипотез в тестах, потому что ящики с усами похожи.
Сейчас я буду применять t-test для сравнения признака INCOME в группах разделения по температуре. У меня сбалансированный дизайн (в группах по 38 городов), поэтому не важно, равны дисперсии или нет, но я проверю. Применять тест Фишера для проверки равенства дисперсий можно только для нормально распределенных (внутри групп) признаков.
hot <- (data_logged %>%
filter(F_TEMPER == "hot"))$INCOME
cold <- (data_logged %>%
filter(F_TEMPER == "cold"))$INCOMEПроверяю нормальность в группах тестом Пирсона (хи квадрат) и Шапира-Уилка
cat(shapiro.test(hot)$p.value, shapiro.test(cold)$p.value)## 0.8588537 0.3617601cat(pearson.test(hot)$p.value, pearson.test(cold)$p.value)## 0.7321958 0.4158802Тесты не обнаружили ненормальности в распределениях в группах. Поэтому применю тест Фишера.
var.test(hot, cold)##
## F test to compare two variances
##
## data: hot and cold
## F = 1.0373, num df = 33, denom df = 31, p-value = 0.9211
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.5101231 2.0928364
## sample estimates:
## ratio of variances
## 1.037312Тест Фишера не отвергает гипотезу о равенстве дисперсий.
Наконец применю t-test, никак не используя результаты про наверно (гипотезу нельзя принимать) равенство дисперсий, потому что дизайн сбалансирован. При нормальных распределениях тест точный.
t.test(cold, hot)##
## Welch Two Sample t-test
##
## data: cold and hot
## t = 0.16315, df = 63.88, p-value = 0.8709
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1722.032 2028.311
## sample estimates:
## mean of x mean of y
## 25361.38 25208.24t-test не отвергает гипотезу о равенстве матожиданий. Применю тест Вилкоксона. Он ранговый и проверят другую гипотезу про сумму рангов.
wilcox.test(cold, hot)##
## Wilcoxon rank sum exact test
##
## data: cold and hot
## W = 569, p-value = 0.7548
## alternative hypothesis: true location shift is not equal to 0Тест Вилкоксона не обнаружил разницы в распределениях, так как не отверг нулевую гипотезу.
В тесте Колмогорова-Смирнова есть смысл, если прошлые тесты не нашли различий. Тест Колмогорова-Смирнова сравнивает функции распределения целиком и проверяет, похожи ли формы распределений.
ks.test(cold, hot)##
## Exact two-sample Kolmogorov-Smirnov test
##
## data: cold and hot
## D = 0.15993, p-value = 0.7266
## alternative hypothesis: two-sidedp-value большое, тест не нашел заметной разницы между формами распределений.
8 сравнение зависимых выборок
Применяю t-test для зависимых выборок.
t.test(
data_logged$POP80,
data_logged$POP92,
paired = TRUE
)##
## Paired t-test
##
## data: data_logged$POP80 and data_logged$POP92
## t = -4.7135, df = 75, p-value = 1.098e-05
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -0.14988791 -0.06083053
## sample estimates:
## mean difference
## -0.1053592Благодаря тому что тест парный, он смог отвергнуть гипотезу о равенстве матожиданий.
Статистика зависимого критерия больше статистики независимого (при положительной корреляции), поэтому мощнее.
Парный тест Вилкоксона
wilcox.test(data_logged$POP80, data_logged$POP92, paired = TRUE)##
## Wilcoxon signed rank test with continuity correction
##
## data: data_logged$POP80 and data_logged$POP92
## V = 730, p-value = 0.0001492
## alternative hypothesis: true location shift is not equal to 0Из-за того что тесты парные, они обнаруживают разницу в распределениях.
III Об анализе зависимостей
1 вспомним pairs plot
if (interactive()) pdf("ggpairs_logged.pdf")
ggpairs(
data_logged[, -(1:2)],
lower = list(continuous = wrap("points", alpha = 0.5, size = 0.3)),
diag = list(continuous = "barDiag")
)
if (interactive()) dev.off()if (interactive()) pdf("ggpairs_logged_colored.pdf")
ggpairs(
data_logged[, -(1:2)],
aes(color = F_TEMPER),
lower = list(continuous = wrap("points", alpha = 0.5, size = 0.6)),
diag = list(continuous = "barDiag"),
columns = names_interesting
) +
scale_color_manual(values = c("blue", "red")) +
scale_fill_manual(values = c("blue", "red"))
if (interactive()) dev.off()2 Начинать нужно с анализа линейных зависимостей
Матрица корреляций пирсона
data_logged <- data_logged %>%
relocate(TEMPER, AREA, POP92, POP80, POPDEN, BORN_F, UNEMP, POVERT, CRIME, INCOME)
if (interactive()) pdf(file = "cor_matrix_pearson.pdf")
ggplot(melt(cor(data_logged %>% select(-CITY, -STATE, -F_TEMPER), method = "pearson", use = "pairwise.complete.obs"))) +
geom_raster(aes(x = Var2, y = Var1, fill = value)) +
geom_text(aes(x = Var2, y = Var1, label = round(value, 2))) +
scale_fill_gradient2() +
theme_dark() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
if (interactive()) dev.off()3 Затем можно переходить к ранговым коэффициентам корреляции.
Матрица корреляций спирмана
if (interactive()) pdf(file = "cor_matrix_spearman.pdf")
ggplot(melt(cor(data_logged %>% select(-CITY, -STATE, -F_TEMPER), method = "spearman", use = "pairwise.complete.obs"))) +
geom_raster(aes(x = Var2, y = Var1, fill = value)) +
geom_text(aes(x = Var2, y = Var1, label = round(value, 2))) +
scale_fill_gradient2() +
theme_dark() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
if (interactive()) dev.off()4 Проинтерпретируйте найденные корреляции - можно ли сказать, что является причиной, что следствием
частные корреляции
Вычитаю из иностранных и плотности доход
(data_logged %>%
dplyr::select(BORN_F, POPDEN, INCOME) %>%
drop_na() %>%
pcor(method = "pearson"))$estimate["BORN_F", "POPDEN"]## [1] 0.5507378(data_logged %>%
dplyr::select(BORN_F, POPDEN, TEMPER) %>%
drop_na() %>%
pcor(method = "pearson"))$estimate["BORN_F", "POPDEN"]## [1] 0.5038904Частная корреляция иностранно рожденных и плотности населения за вычетом дохода больше чем обычная корреляция.
if (interactive()) pdf(file = "pcor.pdf")
ggplot(melt((data_logged %>%
dplyr::select(UNEMP, POVERT, INCOME, CRIME) %>%
drop_na() %>%
pcor(method = "pearson"))$estimate), aes(x = Var2, y = Var1)) +
geom_raster(aes(fill = value)) +
geom_text(aes(label = round(value, 2))) +
scale_fill_gradient2()
if (interactive()) dev.off()
if (interactive()) pdf(file = "cor.pdf")
ggplot(melt(cor(data_logged %>%
dplyr::select(UNEMP, POVERT, INCOME, CRIME) %>%
drop_na())), aes(x = Var2, y = Var1)) +
geom_raster(aes(fill = value)) +
geom_text(aes(label = round(value, 2))) +
scale_fill_gradient2()
if (interactive()) dev.off()
if (interactive()) pdf(file = "pairs.pdf")
ggpairs(
data_logged %>%
dplyr::select(UNEMP, POVERT, INCOME, CRIME) %>%
drop_na(),
lower = list(continuous = wrap("points", alpha = 0.5, size = 0.3)),
diag = list(continuous = "barDiag")
)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
if (interactive()) dev.off()